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Cell colonies of bacteria, tumour cells and fungi, under nutrient limited growth conditions, exhibit complex 
branched growth patterns. In order to investigate this phenomenon we present a simple hybrid cellular automa- 
ton model of cell colony growth. In the model the growth of the colony is limited by a nutrient that is consumed 
by the cells and which inhibits cell division if it falls below a certain threshold. Using this model we have inves- 
tigated how the nutrient consumption rate of the cells affects the growth dynamics of the colony. We found that 
for low consumption rates the colony takes on a Eden-like morphology, while for higher consumption rates the 
morphology of the colony is branched with a fractal geometry. These findings are in agreement with previous 
results, but the simplicity of the model presented here allows for a linear stability analysis of the system. By 
observing that the local growth of the colony is proportional to the flux of the nutrient we derive an approximate 
dispersion relation for the growth of the colony interface. This dispersion relation shows that the stability of 
the growth depends on how far the nutrient penetrates into the colony. For low nutrient consumption rates the 
penetration distance is large, which stabilises the growth, while for high consumption rates the penetration dis- 
tance is small, which leads to unstable branched growth. When the penetration distance vanishes the dispersion 
relation is reduced to the one describing Laplacian growth without ultra-violet regularisation. The dispersion 
relation was verified by measuring how the average branch width depends on the consumption rate of the cells 
and shows good agreement between theory and simulations. 

PACS numbers: 87.18.Hf, 61.43.Hv 



I. INTRODUCTION 

Pattern formation in living systems has attracted much at- 
tention since the pioneering work of D'Arcy Thompson (]]]. 
In recent years special attention has been given to patterns 
emerging from cell colony growth in hostile environments 
llH U . These systems tend to exhibit complex growth pat- 
terns when the growth is limited by the diffusion of a nutrient 
that is necessary for the growth of the cells. The morpholo- 
gies obtained from these living systems resemble that of many 
non-living systems like electrodeposition [4], crystal growth 
dit] and viscous fingers 10, 01- All of these non-living sys- 
tems obey the same underlying growth principle, which is that 
of Laplacian growth, in which the interface between the two 
phases is advanced at a rate proportional to the gradient of a 
potential field. In the case of electrodeposition it is the electric 
field around the substrate, in crystal growth it is the tempera- 
ture field and in viscous fingering the pressure in the liquid. 
This similarity between biological and non-living diffusion 
limited patterns has led to the hypothesis that the biological 
patterns could be explained with the same basic principles [8]. 

Perhaps the most studied example of cell colony growth is 
the growth of bacterial colonies subject to low nutrient levels. 
Bacteria are usually grown in petri dishes at high nutrient con- 
centration. These conditions give rise to colonies with simple 
compact morphologies, but when the growth occurs in more 
hostile low nutrient concentrations the morphologies of the 
colonies take on very complex shapes. This phenomena was 
first reported by Matsushita et al. [8] and since then several 
models have been suggested to explain the observed mopholo- 
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gies. The main modelling approach that has been used is to 
consider the grow th via a system of reaction-diffusion equa- 
tions HI [T(| 1 1 lL [l2ll . These models are able to reproduce 
the observed patterns, ranging from Eden-like II 1311 and dense 
branched morphologies lfl4ll to DLA-like patterns lfl5ll . An- 
other approach by Ben-Jacob et al. iflrj Wfa is to model the 
bacteria as clusters of discrete walkers which obey dynami- 
cal rules. This model also agrees well with experimental data 
and is more biologically realistic compared to the reaction- 
diffusion approach. 

Avascular tumours also grow under similar nutrient limited 
conditions as bacteria cultured in petri dishes. In the early 
stages of cancer development the tumour has yet to acquire its 
own vasculature and the cancer cells therefore have to rely on 
diffusion as the only means of nutrient transport [18]. When 
the tumour reaches a critical size the diffusion of nutrients 
is not enough to supply the inner parts of the tumour with 
oxygen, this leads to cell death or necrosis in the core of the 
tumour. Surrounding the necrotic core is a rim of quiescent 
cells and on the outer boundary there is a thin rim of prolifer- 
ating cells. The mitotic activity therefore only takes place in a 
small fraction of the tumour, while the majority of the tumour 
consists of cells that are either quiescent or dead. Although 
the growth of a tumour is a much more complex process com- 
pared to the growth of bacteria in petri dishes there is still 
evidence from both ex peri ments 11191 |2(J, 12111 and mathemati- 
cal models i22l l23l l24l l25ll that tumours exhibit fingering and 
fractal morphologies driven by diffusion limited growth. 

Another biological system that displays complex patterns 
under diffusion limited growth are fungal colonies. Com- 
plex patterns with fractal morphologies have been observed 
for both multi-cellular filamentous growth 12611 and for yeast- 
like unicellular growth 12711 . These patterns arise in low nu- 
trient conditions or when there is a build up of metabolites 
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which inhibit the growth of the colony and have successfully 
been modelled using both continuous li28l l29tl and discrete 
EH techniques. 

Bacterial colonies exhibit branches which have a width of 
approximately 0.5 mm, which is of the order of 100 cells iflltl . 
This is very different from viscous fingers for example, where 
the disparity of length scales between the molecules and pat- 
tern is much larger. We believe that in order to fully capture 
the dynamics of such systems, where the characteristic length 
scale of the pattern is similar to that of the cells which consti- 
tute the pattern, it is necessary to model them at the level of 
single cells. In this paper we therefore present a simple hybrid 
cellular automaton model of nutrient dependent cell colony 
growth where each automaton element represents a single cell. 
The aim of this model is not to represent any specific biolog- 
ical system, but rather to show that complex growth patterns 
can emerge from a very simple model with minimal assump- 
tions about the cell behaviour. The simplicity of the model 
presented in this paper ensures both generality of the results 
discussed as well as allowing us to carry out a stability analy- 
sis. This analysis we hope will shed light on the growth insta- 
bilities observed in the above mentioned systems. 



II. THE MODEL 

The domain of the colony is restricted to a 2-dimensional 
substrate and the growth is assumed to be limited by some 
essential nutrient which is required for cell division. The sub- 
strate on which the cells grow is represented by aN x N cel- 
lular automaton with lattice constant A. Each automaton ele- 
ment can be in three different states: (i) empty, (ii) hold an ac- 
tive cell or (iii) hold an inactive cell and each element is iden- 
tified by a coordinate x — A(i,j) i,j = 0,l,2,...,N—l. The 
cellular automaton is coupled with a continuous field c(x,t) 
that describes the nutrient concentration. In the case of bac- 
teria the nutrient represents peptone, for tumours oxygen and 
for fungi some sort of carbon source like glucose. The tran- 
sition from an active cell to an inactive occurs if the nutrient 
concentration falls below some critical threshold c„. This in- 
active state corresponds to sporulation or cell quiescence and 
is assumed to be irreversible. An active cell divides when 
it has reached maturation age, it then places a daughter cell 
at random in an empty neighbouring grid point (using a von 
Neumann neighbourhood). If none of the neighbouring grid 
points are empty the cell division fails, but the cell remains 
in the active state. After cell division has occurred the age of 
the parent cell is reset, which means that is has to reach mat- 
uration age again to divide. In order to account for variation 
in maturation age between different cells the maturation age 
of each cell is chosen randomly from a N(z,a) normal dis- 
tribution, where x represents the average maturation age and 
the variance is set to o = x/2. For simplicity we will con- 
sider non-motile cells (which for bacteria corresponds to high 
agar concentrations 081]), which implies that the growth of the 
colony is driven by cell division. 

Active cells are assumed to consume nutrients at some fixed 
rate k, while inactive cells don't consume any nutrients. The 



nutrient is assumed to diffuse in the substrate with a diffusion 
constant D. The nutrient concentration field therefore obeys 
the equation, 



dc(x,t) 



— DV 2 c(x,f) - kn(x,t) 



(1) 



where n(x,t) = 1 if the automaton element at x holds an ac- 
tive cell and n(x,t) — if it holds an inactive cell or is empty. 
The boundary conditions satisfied by the nutrient fields are 
Dirichlet with a constant value c„. This represents a contin- 
uous and fixed supply of nutrient from the boundary of the 
system. In order to simplify the system we introduce new 
non-dimensional variables given by, 
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Using these new variables the equation describing the nutri- 
ent concentration becomes (omitting the primes for notational 
convenience) 



dc(x,t) 



V 2 c(x,t) —kn(x,t). 



(3) 



This equation is discretised using standard five-point finite 
central difference formulas and solved on a grid with the same 
spatial step size as the cellular automaton. Each time step of 
the simulation the nutrient concentration is solved using the 
discretised equation and all the active cells on the grid are up- 
dated in a random order. 



A. Simulations 

Using this simple model of cell colony growth we have in- 
vestigated how the nutrient consumption rate k of the cells 
affects the growth dynamics of the colony. Note that varying 
the non-dimensional consumption rate k is equivalent to ei- 
ther varying the dimensional consumption rate or the bound- 
ary concentration Co», see eq. (f2]). All simulations were started 
with an initial circular colony (with radius 10 cells) of active 
cells at the centre of the grid and an initial homogeneous nu- 
trient concentration of c(x,t) = 1. Figure Q] shows the colony 
after 300 cell generations for k = 0.01,0.1 and 0.2 with x = 10, 
c„ = 0.1 on a grid of size N = 800. From this figure it is obvi- 
ous that the consumption rate of the cells affects the morphol- 
ogy of the colony. For the lowest consumption rate k = 0.01 
the colony grows with a compact Eden-like morphology. The 
colony consists mostly of inactive cells with an outer rim of 
active cells at the boundary. For k = 0.1 the morphology is 
no longer compact but exhibits a pattern similar to the dense 
branching morphology (DBM) observed in viscous fingering 
0. Again the colony consists mostly of inactive cells and 
the few active cells reside on the tips of the branches. For the 
highest consumption rate k = 0.2 the branched morphology 
is even more evident and it exhibits thinner branches. In or- 
der to characterise the morphologies further we measured the 
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fractal dimension of the colonies by measuring how the num- 
ber of cells N (active and inactive) depends on the radius R of 
the pattern Bill . For a compact morphology we expect that 
N ~ R 2 , which is what we find for k = 0.01, but for the two 
other morphologies we find that N ~ Br, where D w 1 .91 for 
k = 0.1 and D « 1 .83 for k = 0.2. For both k = 0. 1 and 0.2 the 
colony thus grows with a fractal morphology. This was also 
confirmed by measuring the density-density correlation func- 
tion C(r) = (p(r')p(r + r')) for the colonies (see inset of fig. 
[T]b and c). At small length scales C(r) decays as r~ a , where 
the fractal dimension of the colony is given by D = 2 — a. For 
k = 0.1 we find a = 0.10 and for k = 0.2 we have a = 0.16, 
which gives fractal dimensions in close agreement with the 
previous method. 

The Eden-like growth pattern observed for k = 0.01 is to 
be expected, as all cells on average divide at uniform speed, 
but what is interesting is that as the nutrient consumption rate 
is increased it leads to a branched morphology. The intu- 
itive explanation of why this type of growth occurs is be- 
cause the high nutrient consumption can't sustain the growth 
of a smooth colony boundary. If a cell on the boundary di- 
vides and places the daughter cell outside the existing colony 
boundary it reduces the chances of neighbouring cells to di- 
vide, as the daughter cell "steals" nutrient at the expense of 
the cells that are closer to the centre of the colony, effectively 
creating a screen from access to the nutrients. It is this screen- 
ing effect that amplifies perturbations to the colony interface 
and leads to the branched morphology. This implies that the 
branched colony morphology is a result of nutrient limited 
growth, which is in agreement with the previously discussed 
experiments and models of colonies of bacteria, tumour cells 
and fungi. 

The dynamics of this model are essentially that of a diffu- 
sion limited growth process, in which the diffusing nutrient is 
transformed by the cells into biomass in the form of daughter 
cells. It is therefore not surprising that it exhibits morpholo- 
gies similar to non-living diffusion limited growth phenomena 
like viscous fingers, electrodeposition and crystal growth. In 
the next section we will quantify the growth instabilities ob- 
served in the system by performing a linear stability analysis 
on the model. 



III. STABILITY ANALYSIS 

In order to analyse the stability of the discrete model we 
have to construct an analogous model that captures the essen- 
tial dynamics but is amenable to mathematical analysis. 



FIG. 1 : Cell colony plots after 300 cell generations for nutrient con- 
sumption rates: (a) k = 0.01, (b) k = 0. 1 and (c) k = 0.2. The other 
parameters were fixed at 1 = 10, c„ =0.1 and grid size N = 800. 
Empty CA elements are coloured white, inactive cells grey and ac- 
tive cells black. This shows that the colony morphology depends on 
the nutrient consumption rate, where a high consumption rate gives 
rise to a fractal colony morphology. The insets in (b) and (c) shows 
a loglog-plot of the density-density correlation function, which show 
that at small length scales C(r) ~ r~ a in the fractal growth regime. 



A. Sharp interface model 

The analogous model will be constructed by considering 
the colony boundary as a sharp interface that moves in two di- 
mensions with a velocity v p determined by the maturation age 
x and the size of the individual cells A = 1, such that v p = 1/x. 
The nutrient consumption of the cells is taken to be k in the 
active part of the colony, where c > c„, zero in the inactive 
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part, where c < c„, and zero outside the colony. If we con- 
sider the growth of a planar front, growing in the y-direction 
and stretching infinitely in the x-direction, the equations de- 
scribing the nutrient concentration are given by, 



dc(x,t) 



V 2 c(x,t), y>y t 



dc(x,t) .. , . . 

I ' = V 2 c(xj)-kH{c(x,t)-c„) 
at 



(4) 



y < yi (5) 



where H is the Heaviside step-function and y, is the position 
of the interface along the y-axis. In order to make the analysis 
simpler we make a change of coordinates to a moving frame 
that travels along with the interface, i.e. we define a new co- 
ordinate £ = y — v p t , where \ = corresponds to y = y,-, the 
position of the interface. This change of coordinates plus the 
fact that there is no dependence on x reduces (@]-[5) to a system 
of three ODE's given by, 



c" + v p c' = 0, E,>0 



c +v„c 



■k=0, -d<t,<0 



;" + V p c' = 0, \<-d 



(6) 



(7) 



(8) 



where ((6]) describes the nutrient concentration outside the 
colony (0 in the active region of the colony and ([8]l in the 
inactive region, and where the width of the active region d 
is determined from the solution of the ODE's. The bound- 
ary condition for the nutrient concentration at \ = °° is that it 
should reach the boundary value c„ — 1 . Moreover we want 
the solution to be smooth across the interface so we require 
that the solutions to (|6]l and (0 have the same value as do 
their derivatives at Z, = 0. We also require that the solutions to 
(01 and ([H) take the value c„ at Z, = —d and that the derivative 
is zero at that point. If we let c e (t,) be the external solution, 
c a (^) the solution in the active region and c/(£) in the inactive 
region we formally require that, 



c e (£) -> 1, as % 
c e (0)=c a (0) 

4(0) = c'M 

Ca(-d) =Cf(-d) 



(9) 



c' a (-d)=c' i (-d)=0. 
A solution to (|6]-[8]) with boundary conditions (O is given by 



k 



^>0 



-d < \ < 



(10) 



where d = (1 — c„)v p /k is the thickness of the boundary layer. 
An example of the solution with appropriate parameter values 
can be seen in fig. [2] where the circles represents the nutrient 
profile measured radially in a simulation with corresponding 
parameter values. The agreement between the two nutrient 
profiles shows that the planar front approach approximates the 
nutrient profile very well. 



B. Interface velocity 

We will shortly analyse the stability of this simple model, 
but before doing so we will discuss the growth dynamics of the 
cells in more detail. In the hybrid model the cells divide at a 
uniform speed, only varied by the stochastic difference in mat- 
uration age x. Although this is the case the model still gives 
rise to interesting growth patterns. The reason behind this is 
that the cells in the model become inactive if the nutrient con- 
centration falls below the threshold c„. If a cell on the bound- 
ary becomes inactive before cell division occurs the interface 
velocity becomes zero at that point, and the inactive cell may 
become the starting point of the development of a fjord. This 
scenario is only possible if the nutrient consumption rate is 
sufficiently high compared to the nutrient concentration at the 
boundary. If this is the case the cells on the boundary have to 
rely on the flux of nutrient in order to survive long enough to 
go through cell division. Our interpretation of this is that in 
the limit of high consumption rates the velocity of the inter- 
face becomes proportional to the flux of nutrient, a mechanism 
already suggested by Matsushita & Fujikawa [8]. Mathemat- 
ically this means that the local interface velocity is given by 
v(x) oc n ■ Vc, where n is the normal of the interface. This 
observation will be the basis for our stability analysis, which 
means that our treatment of the system will not be rigourously 
related to our model, but rather aimed more at understanding 
the dynamics of the model from a qualitative point of view. 



C. Instability of the interface 

In the above solution dTOb of (|6] - © we assumed that the 
interface at \ = was flat, we now introduce a small os- 
cillating perturbation of amplitude 8 to the interface, giving 



C(Q 0.6 - 




FIG. 2: The nutrient profile plotted in the moving ^-frame for k = 
0.03, c n = 0.1 and T = 10 in (?) the inactive part (if) the active part 
and (Hi) outside the colony. The circles represent the radial nutrient 
profile from a simulation and the solid line is the analytic nutrient 
profile (Toll. 
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t,o(x,t) — S(f)cosgx, where 8 C 1. This changes the nutri- 
ent field in the vicinity of the interface, and we need to find 
this perturbed field cg(x,^) to determine the stability. In the 
following analysis we will assume that the interface veloc- 
ity v p <C 1, which means that there is a separation in the 
time-scales between the movement of the interface and the 
dynamics of the nutrient field. This allows for a number of 
simplifications: Firstly it implies that the nutrient field is in 
a quasi-stationary state, which means that the nutrient con- 
centration approximately satisfies V 2 c = outside the colony 
and implies that we can approximate the nutrient profile by 
a linear function in the vicinity of the interface. This is gen- 
erally the case for the types of biological system discussed 
here, where the nutrient diffuses on a time-scale much faster 
than the growth of the colony. For example the reproduction 
time of bacterial cell is of the order of hours, while the dif- 
fusion constant of nutrient in agar is of the order 10 -7 cm 2 /s 
ll32ll . This means that the diffusion time across one bacteria is 
Af « 0.1 s (assuming that a bacteria is approx. 10 /jm), which 
is considerably smaller than the reproduction time. Cancer 
cells are approximately 25 /jm in diameter and the diffusion 
constant of oxygen in tissue is 1 .8 x 10~ 5 cm 2 /s ll33ll . giving a 
diffusion time of Af » 4 x 10~ 3 s across one cell, which is sev- 
eral orders of magnitude smaller than the reproduction time of 
a cancer cell, which is of the order of 10-20 hours. 
Secondly the quasi-stationary assumption allows us to omit 
any time dependence in the solutions for the perturbed field. 
It also implies that the iso-concentration curve — d{x) defined 
by c(x,cj) = c„ will be stationary in the moving frame. Fur- 
ther we will assume that this curve is given by displacing the 
interface by — d in the ^-direction, i.e d(x) =d — Scosgx (cf. 
fig. |3). This is of course only valid when d is small and when 
the wave number q of the oscillation is small. The values of 
d which give rise to branching patterns are of the order of one 
cell size and the interesting range of wave numbers will be 
small as we are not interested in perturbations of wave length 
smaller than a cell size (q < 2k). This means that this assump- 
tion will be valid within the dynamically interesting range. 

The equation of the perturbed nutrient field can now be 
written as 



where the linear part c(^) is given by 



SBe-^ +d) cosqx. 



(11) 



6{Q = l-k/v p (l-e-^)+Wl-k/vj(l-e-^)-c n ) 

(12) 

and B = c'(^) = constant is determined from the boundary 
condition c§(x, — d{x)) — c„. This field satisfies V 2 c = and 
the boundary condition cg(x, — d(x)) = c n (to first order in 8) 
and is therefore an approximate solution for the perturbed in- 
terface. 

As the nutrient field now depends on x the growth of the 
interface is as argued above proportional to n ■ Vcg(x,cjo(x)), 
where n = (1 + S 2 g 2 sin 2 gx)~ 1//2 (8gsingx, 1). But as 8 <C 1 
the interface velocity in the x-direction is negligible and the 
gradient dependent growth velocity can be approximated by 



v{x) 



= A (c , (^(x)) + 8B^-^ Sco " / - ,c+ ' /) cos ? x) (13) 

where A is the constant of proportionality. The velocity can 
also be written as 



v(x) = = 8(f) cos gx. 
at 



(14) 



Taking the derivative in the x-direction and equating the two 
expressions for the velocity gives (only taking into account 
first order in 8) 



d 2 ^ d 2 



C'6 



dtdx dtjdx 



5=5oW 



-Sqsinqx = -AB5q 2 e- q ^°° sqx+ ^ siaqx » -AB8q 2 e- qd sinqx. 
The growth rate 8/8 of the perturbation is therefore given by 
©($)= 8/5 =ABqe- qd ~ 



> kdqe 



-dq 



(15) 



From this we can see that the thickness of the boundary layer 
d affects the stability of the interface. The wave number which 
has the highest growth rate is 



Qmax — 1 1 & 



(16) 



and when d is large (k is small) only modes with a small 
q (long wavelength) have a significant growth rate, but for 
smaller d (larger k) the maxima is shifted to larger q (smaller 
wavelengths) and the growth rates of these wavelengths in- 
crease (cf. fig. m. 

Qualitatively the dispersion relation ( fTBI l can be explained 
in the following way: A perturbation of the colony interface 
gives rise to an identical perturbation in the isoconcentration 
curve c = c„. As the perturbed field is quasi-stationary the 
perturbations decay exponentially in the direction of the in- 
terface ( TTTb . The larger the distance d is between the curve 
c = c„ and the interface, the more the perturbations in the nu- 
trient field decay. Since the interface velocity is proportional 
to the gradient of the nutrient field this implies that the thicker 
the boundary layer is the more uniform the interface velocity 



%(x) = S cos qx 



c s (x£) 




5f=5oW 



FIG. 3: This figure shows the structure of the interface. It is assumed 
that the curve cg(x,^) = c„ is given by displacing the inteface by — d 
in the ^-direction. 
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is, and consequently the interface is less sensitive to perturba- 
tions. 

In the case d = the dispersion relation reduces to ($(q) ~ 
q, which is the dispersion relation for Laplacian growth with- 
out ultra-violet regularisation 113411 - In this type of growth the 
interface velocity is proportional to a potential field which is 
zero on the interface. In the above derivation the field which 
governs the growth of the interface instead takes on a zero 
value at a distance d from the interface. With the dispersion 
relation of Laplacian growth in mind we can conclude that 
by introducing a boundary layer the interface is stabilised for 
high wave numbers, but that this stabilising effect decreases 
as the thickness of the boundary layer is reduced. As men- 
tioned before the thickness of the boundary layer d ~ 1/k, 
which means that the stability of the colony growth depends 
directly on the consumption rate of the individual cells. A low 
nutrient consumption results in a wide boundary layer that sta- 
bilises the growth, while a high consumption gives rise to a 
thin boundary layer which leads to unstable branched growth. 

The reason why all wave numbers have a positive growth 
rate ((o(q) > for all q) is because the dynamics do not con- 
tain any stabilising mechanism. In the Mullins-Sekerka insta- 
bility [35], which also describes a diffusion-limited growth, 
the growth is stabilised by surface tension, which inhibits the 
growth of perturbations with a large wave number. A similar 
type of effect can be observed in a reaction-diffusion model 
describing bacterial growth Il36ll . In the reaction-diffusion 
model a protrusion into the nutrient side of the interface re- 
sults in enhanced local growth, but the bacterial diffusion flow 
is is reduced at the protrusion. It is the relative strength be- 
tween these two effects that determines the stability of the 
growth. This type of stabilisation does not occur in our sys- 
tem because the cells are immobile and the growth does not 
depend on the local curvature of the interface. Consequently 
there are no perturbations that have negative growth rate. 

Another system which exhibits a Mullins-Sekerka like in- 
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FIG. 4: This plot shows the dispersion relation j 1 5b for k = 
0.03,0.1,0.2 and x= 10, c„ =0.1. It can be seen that both the fastest 
growing mode and its growth rate depends on the consumption rate 
k. 



stability is the Fisher equation with a cut-off in the reaction 
term for low concentrations of bacteria l37ll . This is mo- 
tivated by the fact that bacteria are discrete entities, which 
means that at some small concentration the continuum formu- 
lation breaks down. Because we consider single cells rather 
than concentrations, the cut-off effect is already incorporated 
in our model. On the other hand we also have a cut-off in the 
growth due to low nutrient concentrations, i.e. no cells divide 
in regions where c < c„. Although this cut-off is due to cells 
becoming inactive rather than finite particle numbers it might 
have a similar effect on the stability of the continuous model 
and is a question certainly worth investigating. 

D. Comparison to Simulations 

The above derivation of the dispersion relation ( fTBI ) con- 
tains a number of simplifications and assumptions and it is 
therefore important to verify the analytic result by compar- 
ing it to simulation results from the discrete model. This was 
done by measuring the average branch width in the colony and 
how it depends on the consumption rate k. The consumption 
rate affects the branching of the colony as it determines the 
linear stability of the interface. When a branch grows it is 
constantly subject to perturbations and when it reaches a crit- 
ical width it becomes linearly unstable and splits, similar to 
what occurs in splitting of viscous fingers 13811 . As we don't 
have any detailed information about the dynamics of the tip 
splitting we will consider a idealised version of the process. 
We will assume that the branches grow to the critical width 
h = \nax = 2.K/q max at which they split and that each split- 
ting gives rise to two branches of equal width. If we assume 
that no branches are annihilated and that they grow at a con- 
stant speed then an estimate of the average branch width in the 
colony is 

lavg ~ + k max )/2 = 3/4\ max = 3/2nd. (17) 

This is of course a highly idealised picture of the branching 
process, but contains the essential dynamics as it is clear from 
figure [TJ that the branch width decreases with increasing k. 

The results from the simulations were analysed in the fol- 
lowing way: First the colony was allowed to grow long 
enough for the morphology to be properly established (ap- 
prox. 400 time steps), the cell density was then measured on 
a circle of radius R = 0.15r nuix , where r max is the distance to 
the most distant cell in the colony (cf. fig. [5]). The resulting 
cell density was stored in a vector n(9), where n(8) = 1 if the 
automaton element at distance R and angle 8 from the centre 
holds a cell (active or inactive) and n(9) = if it is empty and 
from this vector the average branch width was calculated. In 
order to make sure that the measurements were not biased by 
the choice of radius we also measured how the average branch 
width depended on the radius. The results show that the aver- 
age branch width depends on the radius for small R, but that 
this bias is negligible for the values of R we consider (data not 
shown). 

The average branch width was calculated for several values 
in the range of k (averaged over 20 simulations for each value 
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of k) where branching occurs and the results can be seen in 
fig. [6] where it is plotted together with the analytic result ( fT71 ). 
From this we can see that the average branch width from the 
simulations agree with the analytic result obtained from the 
linear stability analysis of the model. The agreement is not 
perfect but the simulation results exhibit an approximate 1 jk 
decay of the average branch width predicted by the stability 
analysis. One should also bear in mind that our analysis con- 
tains a number of simplifications which means that we cannot 
expect to capture the exact dynamics of the system, but at least 
our analysis predicts the general behaviour of the system. 

IV. CONCLUSION 

In this paper we have presented a simple hybrid cellular 
automaton model of cell colony growth, which exhibits inter- 
esting growth patterns. We have investigated how the nutrient 
consumption rate of the cells (or equivalently the nutrient con- 
centration) affects the growth dynamics. The results show that 
for low consumption rates the colony grows with a Eden-like 
morphology, while for higher consumption rates the colony 




FIG. 5: The circle around which the average branch width was mea- 
sured. 



avg 
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k 

FIG. 6: The average branch width as a function of the consumption 
rate k. The circles shows the average result from simulations and 
the solid line represents the analytic result. The error bars show the 
standard deviation of the simulation results. 



breaks up into a branched morphology. By observing that 
the local growth of the colony is proportional to the gradi- 
ent of the nutrient field we were able to derive a dispersion 
relation, which shows that the thickness of the boundary layer 
in the colony determines the stability of the growth. When the 
nutrient consumption rate is low the colony exhibits a wide 
boundary layer, which stabilises the growth, while when the 
consumption is high the width of the boundary layer is re- 
duced and the growth becomes unstable leading to a branched 
morphology. When the boundary layer vanishes the derived 
dispersion relation is reduced to the one describing Lapla- 
cian growth without ultra-violet regularisation. Analysis of 
colonies obtained from the discrete model show good agree- 
ment between simulations and theory. 

Some cells are known to use chemotactic signalling under 
harsh growth conditions. This has for example been observed 
in bacterial colonies, which under very low nutrient conditions 
exhibit densely packed radial branches [16]. It has been sug- 
gested that this occurs because stressed cells in the interior of 
the colony secrete a signal which repels other cells. This could 
be included in the model either by introducing a bias towards 
placing the daughter cell in the neighbouring automaton ele- 
ment that has the lowest (or highest) level of the chemotactic 
substance or by allowing cells to move down (or up) gradi- 
ents of the substance. Introducing a chemorepellant secreted 
by cells exposed to low nutrient concentrations would most 
likely lead to a more directed growth away from the colony 
centre, and thus to a more ordered morphology with straighter 
branches. Another approach could be to introduce a direct 
chemotactic response to the nutrient, which probably would 
have a similar effect on the colony morphology. It should be 
noted that the introduction of chemotaxis would make the dy- 
namics of the model more complex and would render the sta- 
bility analysis much more difficult. 

The instabilities described in this paper should be observ- 
able in any system where cell colony growth is limited by 
some nutrient that diffuses and which leads to inhibition if 
it falls below some critical threshold. The nutrient field also 
has to be in a quasi-stationary state, which corresponds to a 
separation in time scales between the cell division and the dy- 
namics of the nutrient field. Additionally we require that the 
colony expansion occurs mainly by cell division at the colony 
boundary and not by movement of individual cells. These 
conditions apply to growth of avascular tumours, bacterial 
colonies grown in high agar concentrations, yeast colonies and 
fungal growth. All of these systems exhibit branched or frac- 
tal morphologies under certain growth conditions and these 
growth patterns may be explained by the dispersion relation 
presented in this paper. 
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